# Test with Covariates 

# We consider three covariates:i) Eastern Europe ii) international security event/context iii) change from being in opposition to being in government.

#Eastern Europe is already in the dataset
# we first measure the international security event

#-----------------------------------------------------------------------------
# RDD with Covariates  
#-----------------------------------------------------------------------------
sink(file = "log_file2.txt", type = "output")

require(tidyverse) 
require(rio)      
require(magrittr)
require(dplyr)
require(stringr)
require(ggplot2)
require(tidyr)
require(MASS)
require(stargazer)
require(coefplot)
require(sjPlot)
require(margins)
require(plm)
require(estimatr)
require(lmtest)
require(lme4)
require(pscl)
require(jtools)
library(broom)     
library(rdrobust)   
library(rddensity)
library(rddtools)
library(huxtable)   
library(rdd)
require(haven)
require(plyr)
require(cowplot)
library(data.table)
require(gridExtra)
library(estimatr)
library(readr)
require(grid)
require(ggpubr)


covdata<-read.csv('covdata.csv')

#creating variable for country FE
covdata$state.f = factor(covdata$countryname)
covdata$state.d = model.matrix(~covdata$state.f+0)


################################################################################

#Table A4

#Position: Military Positive - Military Negative Per104-Per105 (with covariates)
covdata$military<-(covdata$per104)-(covdata$per105)
covdata$military_l<-log(covdata$military+1)
covdata_1<-covdata%>%dplyr::select(dif_l1, military_l, East, UCDP_intervention, cabinet_party,countryname,state.d, treatment)
covdata_1<-covdata_1[complete.cases(covdata_1), ]
fakefuzzy <- rdd_data(x=covdata_1$dif_l1, y=covdata_1$military_l,z=covdata_1$treatment,covar=cbind(covdata_1$East,covdata_1$UCDP_intervention,covdata_1$cabinet_party, covdata_1$state.d),cutpoint=0)
fakefuzzy2 <- rdd_data(x=covdata_1$dif_l1, y=covdata_1$military_l,z=covdata_1$treatment,covar=cbind(covdata_1$East,covdata_1$UCDP_intervention,covdata_1$cabinet_party),cutpoint=0)


#parametric 
summary(f1<-rdd_reg_lm(rdd_object=fakefuzzy,covar.opt = list(strategy = c("include"), slope = c( "separate")),order=1))
coeftest(f1, vcov=vcovHC(f1,type="HC0",cluster="countryname"))
summary(f2<-rdd_reg_lm(rdd_object=fakefuzzy2,covar.opt = list(strategy = c("include"), slope = c( "separate")),order=1))
coeftest(f2, vcov=vcovHC(f2,type="HC0",cluster="countryname"))

#non-parametric 
summary(rdbwselect_2014(y = covdata_1$military_l, x = covdata_1$dif_l1, c=0, bwselect="CCT"))
#FE
summary(rdrobust(y = covdata_1$military_l, x = covdata_1$dif_l1, c=0, kernel = "tri",level = 95, p=2, fuzzy = covdata_1$treatment, h=1.844, b=4.057,covs=cbind(covdata_1$East,covdata_1$UCDP_intervention,covdata_1$cabinet_party, covdata_1$state.d),cluster=covdata_1$countryname, all=TRUE))
#w/o FE
summary(rdrobust(y = covdata_1$military_l, x = covdata_1$dif_l1, c=0, kernel = "tri",level = 95, p=2, fuzzy = covdata_1$treatment, h=1.844, b=4.057,covs=cbind(covdata_1$East,covdata_1$UCDP_intervention,covdata_1$cabinet_party),cluster=covdata_1$countryname, all=TRUE))

#Figure A3
d1<-covdata_1%>%dplyr::filter(dif_l1>-5, dif_l1<5)
pdf("Fig.A3.pdf") 
rdplot(x=d1$dif_l1, y=d1$military_l, c=0, p=2,covs=cbind(d1$East,d1$UCDP_intervention,d1$cabinet_party, d1$state.d), x.label="Difference between RRPP vote and threshold (t-1)", y.label="Defence Position" , title="")
dev.off()

sink()
